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Observations of the circumstellar disk in the Bok globule CB 26 at 110, 230, and 270 GHz are 
presented together with the results of the simulations and estimates of the disk parameters. These 
observations were obtained using the SMA, IRAM Plateau de Bure, and OVRO interferometers. 
The maps have relatively high angular resolutions (0.4-1"), making it possible to study the spatial 
structure of the gas-dust disk. The disk parameters are reconstructed via a quantitative comparison 
of observational and theoretical intensity maps. The disk model used to construct the theoretical 
maps is based on the assumption of hydrostatic and radiative equilibrium in the vertical direction, 
while the radial surface density profile is described phenomenologically. The system of equations for 
the transfer of the infrared and ultraviolet radiation is solved in the vertical direction, in order to 
compute the thermal structure of the disk. The disk best-fit parameters are derived for each map 
and all the maps simultaneously, using a conjugate gradient method. The degrees of degeneracy of 
the parameters describing the thermal structure and density distribution of the disk are analyzed 
in detail. All three maps indicate the presence of an inner dust-free region with a radius of approxi- 
mately 35 AU, in agreement with the conclusions of other studies. The inclination of the disk is 78°, 
which is smaller than the value adopted in our earlier study of rotating molecular outflows from CB 
26. The model does not provide any evidence for the growth of dust particles above a max ~ 0.02 cm. 



1. INTRODUCTION 

The rapidly growing number of detected planetary 
systems, which have very different orbital configura- 
tions and orbit stars with masses varying over a wide 
range 1 , suggests that the formation of planets is a 
widespread process in the Galaxy. Since gas-dust 
disks are natural predecessors of planetary systems, 
theoretical and observational studies of such disks 
have become intense over the past decade. 

According to current understanding, stars are 
formed during the gravitational collapse of molec- 
ular clouds Since clouds always have non-zero 
angular momenta, the collapsing matter cannot di- 
rectly fall into a protostar, and fairly rapidly (over 
~ 10 4 yr) forms a circumstellar disk surrounded by 
an envelope |2jj. The angular momentum inside the 
disk is redistributed so that the bulk of the matter 
falls onto the protostar, and only a smaller portion 
moves outward, carrying away angular momentum 
The characteristic evolutionary timescales of 
dust disks around young single stars are several mil- 
lion years [6j. At this time, a key mechanism is ini- 
tiated that can eventually lead to the formation of 
planets: the growth of the size of the dust parti- 
cles and their settling toward the central plane of 
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the disk. The enlargement of dust particles with the 
subsequent formation of planetesimals form the ba- 
sis of the core accretion theory Q , which appears to 
describe the main regime of planetary formation Q . 
The evolution of dust in disks around young stars is 
confirmed by their observed spectral energy distri- 
butions (SEDs) in the IR and (sub) millimeter (9llicj: 
together with the high rate of occurence of exo- 
planets, this suggests a protoplanetary nature for 
these disks. Additional important factors affecting 
the evolution of young stars are bipolar jets, out- 
flows, and disk winds [11], which play an important 
role in the final stages of evaporation of the gaseous 
disk. The variety of controlling processes and the 
wide range of their physical parameters makes proto- 
planetary disks very interesting objects; a complete 
understanding of their nature will be only possible 
through a synthesis of modern observations and de- 
tailed simulations. 

Observations of protoplanetary disks are difficult, 
since these objects have relatively small sizes and low 
temperatures, but a number of basic facts concern- 
ing their evolution and structures have been more 
or less reliably established [T3 |, Observations in 
the middle-IR can be used to determine the rate 
of occurence and lifetime of these disks [l3[ . Their 
masses can be determined using millimeter observa- 
tions [14|], and the structures of protoplanetary disks 
can be reconstructed using interferometric observa- 
tions [15]. A number of tools have been developed 
as a theoretical basis for understanding the physics 



2 



of protoplanetary disks, and for the interpretation 
of SEDs (see, e.g. , the list of publicly available soft- 
ware codes in |16l|). However, the problem of degen- 
eracy remains, when observations are equally well 
described by several sets of parameters, and a model 
can only provide limits for these sets [17| . Increasing 
the angular resolution of observations can help to re- 
solve this problem. Spatially resolved observations 
are currently available for more than one hundred 
and fifty protoplanetary disks at wavelengths from 
the optical to the radio, in both molecular lines and 
the continuum 2 . These high-quality observational 
data require adequate interpretations and new ap- 
proaches to deriving as much information as pos- 
sible about the physical and chemical structure of 
the disk. The successful determination of the pa- 
rameters of a number of protoplanetary disks us- 
ing multi-frequency, spatially resolved observations 
[let [T^ | and the operation of the ALMA interferom- 
eter make this area promising. 

The aim of our study is to determine a self- 
consistent physical structure for the protoplanetary 
disk located at the edge of the Bok globule CB 26. 
This disk is about a hundred AU in size, has a 
mass of O.IMq, and belongs to Class I of young 
stellar objects according to the classification of (20j . 
The disk in CB 26 has been well studied and mod- 
eled [U • A rotating outflow responsible for the 
removal of ang ular momentum from the disk was 
detected in [23j . A central region 45 ± 5 AU in size 
where dust is depleted was later discovered in the 
disk using spatially resolved observations at 230 GHz 
[24| . The best-fit model was derived by compar- 
ing observed images and synthetic maps. However, 
due to the computationally intensive method used 
for the radiative-transfer calculations, a step-by-step 
method was used to determine the parameters, and 
possible degeneracy of the model parameters was not 
analyzed. Another shortcoming of the above work is 
the phenomenological law used for the disk density 
distribution in the vertical direction. Our aim is to 
develop new tools to reconstruct disk structure and 
advanced routine to identify best-fit model param- 
eters and their degeneracies. We aim to determine 
to what degree the modeling of mm maps of CB 26 
are reproducible and reliable in the light of (a) new 
observations with the IRAM Plateau de Bure Inter- 
ferometer (PdBI), (b) a refined model for the vertical 
density distribution, and (c) a more well-grounded 
parameter search technique. 

The paper contains four sections, describing the 
observations, a physical model for the disk, the pro- 
cedures used to obtain the synthetic maps, and our 
searches for the best-fit model. The main results and 
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conclusions are given at the end of the paper. 



2. OBSERVATIONS AND DATA REDUCTION 



2.1. OVRO observations 

CB 26 was observed with the Owens Valley Radio 
Observatory (OVRO) interferometer between Jan- 
uary 2000 and December 2001. The 1.3 and 2.7 mm 
continuum emission was observed simultaneously in 
2-GHz bands using an analog correlator, except in 
the highest-resolution configuration, where the 4- 
GHz capability of the new 1 mm receivers was used. 
Four configurations of the six 10.4 m antennas pro- 
vided baselines in the range 10-170kA at 2.7 mm 
(110 GHz) and 10-370 kA at 1.3 mm (236 GHz). 
The average system noise temperatures of the He- 
cooled SIS receivers were 300-400 K at 110 GHz 
and 300-600 K at 236 GHz. The observing param- 
eters are summarized in Table [T] 

The amplitude and phase calibration were based 
on frequent observations of a nearby quasar. The 
flux densities were calibrated using observations of 
Uranus and Neptune, yielding relative uncertainties 
of 20%. The raw data were calibrated and edited 
using the MM A software package [2^. The map- 
ping and data analysis were carried out using the 
MIRI AD package [26j . The final maps were obtained 
by combining the OVRO and IRAM PdBI data. 

2.2. IRAM PdBI observations 

Observations of CB 26 were carried out with the 
IRAM Plateau de Bure Interferometer in Novem- 
ber 2005 (D configuration with five antennas) and 
December 2005 (C configuration with six antennas; 
project PD0D). Two receivers were used simulta- 
neously, tuned to single side bands at 89.2 GHz 
and 230.5 GHz, respectively. Higher resolution ob- 
servations at 230 GHz were carried out in Jan- 
uary 2009 (B configuration with six antennas) and 
February 2009 (A configuration with six antennas; 
Project S078). Further observations at 86.7 GHz 
were obtained in November 2008 (C configuration) 
and March 2009 (D configuration; project SC1C) 
with six antennas, using the new 3 mm receivers. 
Several nearby phase calibrators were observed dur- 
ing each track to determine the time-dependent com- 
plex antenna gains. The correlator bandpass was 
calibrated using 3C 454.3 and 3C 273, and the ab- 
solute flux-density scale was derived from observa- 
tions of MWC 349. The flux calibration uncertainty 
is estimated to be <20% at both wavelengths. The 
observing parameters are summarized in Table [TJ 

The raw data were calibrated and imaged using 
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Table 1. Millimeter interferometric observations of CB 26 in the continuum: Ao and vq are the central wavelength and 
frequency, 9 the width of the antenna main beam, a and b the semi-major and semi-minor axes of the interferometer beam, 
PA the position angle of the beam major axis measured from the northern protoplanetary disk axis in the clockwise direction, 
a v the rms map uncertainty in mJy/£2„, and Q„ = 1.133a6 the effective solid angle of the beam. 
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a) Composite maps have been rotated by 32° in order for the plane of the disk to align with the x axis. 



the latest version of the GILDAS 3 software. The 
final maps were obtained by combining the OVRO 
and PdBI data. 

2.3. SMA Observations 

Observations with the Submillimeter Array (SMA, 
[22}) were made on December 6, 2006 (extended con- 
figuration) and December 31, 2006 (compact config- 
uration), covering frequencies from 267 to 277 GHz 
in the lower and upper sidebands, respectively, and 
providing baselines of 12-62kA. The typical sys- 
tem temperatures were 350-500 K. The quasar 
3C 279 was used for the bandpass calibration, and 
the quasars B0355+508 and 3C 111 for the gain 
calibration. Uranus was used for the absolute flux 
calibration, which is accurate to 20% -30%. The 
1.1 mm continuum map was constructed using line- 
free channels in both sidebands. The main observa- 
tional parameters are summarized in Table [TJ The 
raw data were calibrated using the IDL MIR pack- 
age [28| and visualized using a package MIRIAD. 



3. PHYSICAL MODEL 

A global aim of this paper is to develop a set 
of software for determining the parameters of pro- 
toplanetary disks using observed millimeter maps. 
The reconstruction of the disk parameters is based 
on a quantitative comparison of synthetic and ob- 
served images of the disk. Mathematically, the 
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problem is reduced to searching for the minimum 
of the parameter-dependent function describing the 
difference between the observed and synthetic maps. 
In this Section, we present the protoplanetary disk 
model used to construct the synthetic maps. 

The main sources of disk heating in the model are 
the radiation of the central star and viscous heat- 
ing (important only in dense central regions), while 
the source of cooling is thermal radiation in the con- 
tinuum. Determining the disk's thermal structure 
is reduced to the problem of radiative transfer in a 
dusty medium, since the dust contributes mostly to 
the opacity of the matter in protoplanetary disks. 
The density distribution in a low-mass disk around 
a single star can be considered to be axially sym- 
metric in a first approximation (this is also valid 
for peripheral regions around close binaries). The 
thickness of the disk increases fairly rapidly with ra- 
dius, so that the central star directly illuminates not 
only the inner edge of the disk, but also its periph- 
eral parts [29]. Moreover, the dust temperature is 
mainly determined by the radiation incident from 
the disk surface, i.e., by the vertical radiation trans- 
fer, since the optical depth rapidly increases in the 
radial direction. This makes it possible to decom- 
pose the two-dimensional problem and construct a 
so-called 1+1D model of the protoplanetary disk, in 
which the vertical and radial structures of the disk 
are determined independently |i30l. l3l|. 

The radial distribution of the surface density in 
our model is given by the analytical expression 

whose parameters are determined from observa- 
tional data. The disk was considered to be hydro- 
static in the vertical direction, and the gas-density 
distribution was calculated based on its tempera- 
ture. The dust and gas were assumed to be well 
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mixed, so that the ratio of their densities did not 
change over the disk, and was equal to 0.01. It was 
additionally assumed that the gas and dust temper- 
atures were equal; this equality is achieved in dense 
regions due to effective collisions between the dust 
particles and molecules. In the disk atmosphere, the 
gas can be much more hotter than the dust [HI], 
but there is only a small mass of dust in this re- 
gion. Therefore, the assumption that the dust and 
gas temperatures are equal is justified when obtain- 
ing synthetic maps of the disk in the mm, where 
the disk is optically thin and the radiation-intensity 
distribution reflects the dust density distribution. 

The calculation of radiation transfer in a proto- 
planetary disk is a computationally intensive task 
due to the high optical depths for photons in the 
UV (up to ~ 10 5 ). Since the computation time re- 
quired to calculate one model is critical to searches 
for the best-fit model, we developed a fast method 
for calculating the dust temperature in the disk with 
acceptable accuracy. This was achieved through an 
accurate transition to a two-frequency approxima- 
tion and the use of the Schuster-Schwarzschild and 
Eddington approximations. We tested the method 
and compared it with another algorithm for calcu- 
lating the radiative transfer, in which the dust opac- 
ity coefficients are functions of the wavelength. The 
transition from monochromatic opacities to mean 
opacities makes it possible to appreciably reduce the 
computational time without losing accuracy when 
calculating the dust temperature. 

3.1. Momentum Equations for the IR Radiation 
Transfer 

Let us consider a ring in the disk at the radius R. 
The gas surface density at this radius, in a layer be- 
tween the equatorial plane and the surface of the 
disk, is Ex- The key assumption of the model is 
that the disk does not radiate in the UV. Thus, we 
divided the spectral range into UV and IR parts and 
describe them separately. It was also assumed that 
the disk is geometrically thin and the plane-parallel 

z 

approximation is valid. If £ — J p(z)dz is the sur- 

o 

face density in a layer from to z, then the momen- 
tum equations for IR radiation transfer in the grey 
and Eddington approximations can be represented 
in the form 
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Here, k„ and K^ ca are the monochromatic coeffi- 
cients of true absorption and scattering. Note that 
Kp and kr are functions of the dust temperature. 

3.2. Balance Equation and the Heating Functions 

The equations considered must be closed by an 
energy-balance equation. The total radiative flux in 
the IR is due to heating of the surrounding medium. 
The change in this flux in the vertical direction is 
described by the formula 
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Here, 5 s tar and S acc [erg s -1 g _1 ] represent the heat- 
ing due to stellar UV radiation and gas accretion. 
The stellar radiation heating function is 



"Sstar — 4:lTKp J uv , 



(7) 



where k p v is the Planck mean true absorption coef- 
ficient in the UV per unit mass, and J uv is the UV 
intensity averaged over angles and frequencies. We 
used the two-flux Schuster-Schwarzschild approxi- 
mation to calculate J uv , assuming that the disk does 
not radiate in the UV, and only absorbs and scatters 
the stellar radiation. The corresponding system of 
equations has the form 



dF uv 



AirKp </ uv , 
47T dJ uv 
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where n F v is the mean flux extinction coefficient (true 
absorption and scattering), and F uv is the UV flux. 

The boundary conditions at the surface of the disk 
can be written in the form 



dF 
F = 



ckp(B-E), 
c dE 



3kr dS ' 



(2) 
(3) 



where F [erg cm~ 2 s _1 ] and E [erg cm -3 ] are the in- 
tegrated energy flux and energy density, B = oT A is 
the integrated density of black-body radiation, and 
Kp and Kp [cm 2 g _1 ] are the Planck and Rosseland 
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(10) 



where J uv is the intensity of the stellar radiation 
averaged over a hemisphere: 



L st ar/(47ri? 2 ). 



(11) 



The coefficient / determines the portion of the stel- 
lar radiation intercepted by the disk. Generally, / 



5 



depends on the inclination of the disk surface relative 
to the incident radiation. It is natural to combine 
the input model parameters / and L stal into a single 
parameter / • £ s tar, representing the portion of the 
stellar luminosity intercepted by the disk. 

The above system of equations can be solved an- 
alytically, if the opacities do not depend on E, or 
numerically, e.g., using the finite difference method. 
The UV opacities are defined as follows: 
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Here, the Planck function B v depends on the stellar 
temperature. Assuming that the gravitational en- 
ergy release is proportional to the density, the heat- 
ing function due to gas accretion is 



Sure- 
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(14) 



where M is the accretion rate, and M sta r is the stel- 
lar mass. 

3.3. Solution of the System of Equations for the 
Thermal Structure 

The above heating functions depend only on the sur- 
face density. Let us assume for the moment that 
Kp and kr are only functions of the surface den- 
sity (they are functions of temperature as well, but 
this problem will be solved further using an iteration 
process). In this case, we can obtain an analytical 
solution for T(E). The temperature can be deter- 
mined from Equations ((2J) and ©, i.e., from the 
relationship 



CKp(aT 4 — E) = S s \ 



S& 



(15) 



It is necessary to find E(E) in order to obtain T(E). 
We integrate Equation © and obtain F = F(E). 
Then, F(E) is substituted in Equation ((3J), and the 

s 

latter is integrated as J dE'. As a result, we obtain 



an equation with the known function <5e(E) 
E(Tt) — E(0) = <5_g(E). 



(16) 



We must know the radiation energy density E(0) in 
the equatorial plane of the disk in order to recon- 
struct E(Y,) from this relationship. The radiation 
density can be found from boundary conditions. If 
the IR radiation from the disk is isotropic at the 



disk surface, i.e., at E=Et, whereas the IR radiation 
from the star can be neglected compared to the ra- 
diation from the disk itself, then -F(Eo) = — oE(Et)- 
Therefore, the radiation energy density at the disk 
surface is 



£(E T ) = 



2F(E T ) 



(17) 



The formulas for E(0) are obtained from ([16]) for 
S = St- 



E{0) = £(E T ) - fe(E T ) 



(18) 



When constructing the solution, we assumed that 
Kp and kr are functions only of the surface den- 
sity, not the temperature. An iterative process was 
implemented to solve this problem. After the ther- 
mal structure was found for the given opacities, we 
recalculated the opacities using information on the 
temperature profile T(E). The new opacities were 
then used to recalculate a refined thermal structure. 
In practice, this iterative process converged over 8 - 
10 steps. 

3.4. Equation of Hydrostatic Equilibrium 

The equation of hydrostatic equilibrium for a geo- 
metrically thin disk has the form 
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where c|ds the sound speed. Since T, and therefore 
c T as a function of E, (but not as function of z) is 
known, the latter equation should be written in the 
form 



dz 
dE 



GAL, 



R :i 



(20) 
(21) 



This system of equations can be solved using the ex- 
plicit integration scheme if the density p(0) in the 
equatorial plane is known; this should correspond to 
the boundary condition p(Ex) = p C xt, where p C xt is 
the given density at the disk surface. To search for 
the required value of p(0), we used a combination 
of the shooting method and bisection method. The 
solution z = z(E) and p = /c(E) of the system of 
Equations (f2"0"|) - (|2"Tj) can be used to obtain the re- 
quired relationships p = p{z) and T — T(z). 

The choice of a grid of spatial coordinates is an 
important part of the model. A logarithmic grid is 
a natural choice for the radial direction. This grid 
smoothly traces the density gradients in the case 
of a power-law parametrization of the surface den- 
sity The left and right boundaries of the radial 
grid are the inner R ln and outer i? ou t radii of the 
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disk. The choice of a grid in the vertical direction is a 
more complex task. If we knew the density distribu- 
tion in the vertical direction a priori, we could place 
the grid nodes so that the density would change by 
no more than some specified factor between adja- 
cent cells. However, the density profile is not known 
a priori. On the other hand, the density distribution 
can be found analytically for a vertically isothermal 
disk with a specified temperature. Therefore, we 
used the grid for an isothermal disk to calculate the 
density distribution in a nonisothermal disk with a 
similar temperature. 

In practice, calculation of disk model on a grid 
of 500x500 cells takes two to three minutes using 
a 3 GHz processor. Figure O shows examples of the 
density and temperature distributions for one model 
of the disk in CB 26. 



4. CONSTRUCTION OF SYNTHETIC MAPS 

After model temperature and density distribu- 
tions in the disk have been constructed, it is neces- 
sary to calculate the spatial distribution of the radia- 
tion intensity I v coming to the Earth for a given disk 
orientation in the plane of the sky, i.e., to construct 
a theoretical image. The theoretical maps were con- 
structed using the tracing routine from the NATALY 
software [33|]. The scattering was neglected, since 
the two-frequency disk model does not enable cal- 
culation of the distribution of the spectral radiation 
intensity J„. Note, however, that scattering is negli- 
gible compared to thermal radiation in the millime- 
ter wavelengths. 

We must know the monochromatic opacities «;f s 
and 4° a to construct theoretical maps and obtain 
the mean opacities in (g}, ©, ©, and (12]). The 
monochromatic opacities were calculated based on 
the true absorption efficiency factor Qabs(^, a) and 
the scattering efficiency factor Q sca (v 7 a), calculated 
using Mie theory: 

a max 

< bB = — [ Q ahs (v,a)ira 2 f(a)da, (22) 
Pd J 

^min 
fl max 

C a = — / Q sc ,(v,a)Tra 2 f(a)da. (23) 
Pd J 

a m i n 

The dust-particle size distribution function f(a) was 
assumed to be a power-law, parametrized in terms of 
the maximum and minimum particle sizes a max and 
a m i n and the slope a pow . An additional parameter 
is the mass fraction of silicate fsi in the mixture of 
silicate and graphite dust particles. 

To obtain a model disk image /* (synthetic map) , 
the ideal disk image I v must be convolved with the 



beam: 



I v (x', y')W(x' - x, y - y) dx'&y' 



(24) 

The full widths at half maximum of the beam 
W(a, b) were determined using the major and minor 
axes HPBW a and HPBW D , and the position angle 
PA: 



W(a,b) 



where 



x exp 
x exp 

HPBW a 



( a cos OL — b sin a 



a sin a-\-b cos a 



(25) 



HPBW fc 



. |_ pa 

2Vln2 2Vln2 2 

(26) 

We focused on the development of a rapid and 
accurate calculation technique, since the convolu- 
tion with the beam is the most resource-intensive 
step of the numerical integration. The multiple inte- 
gral (|2~4"|) reduces to an iterated integral that was cal- 
culated using a composite quadrature formula with 
a highest degree of accuracy (the composite Gaus- 
sian quadrature); i.e., the entire integration interval 
was divided into unequal subintervals, within each 
of which a simple Gaussian quadrature was used. 
The subintervals are searched for using a local dou- 
ble recalculation method: a subinterval is divided 
into two new subintervals only if it yields a signifi- 
cant improvement in the integration accuracy. This 
adaptive control of the integration accuracy enables 
us to appreciably reduce the time required to com- 
pute (j2"4]h 



5. SEARCH FOR THE BEST-FIT MODEL 

We used an analog of the reduced x^ ed as a crite- 
rion for agreement between the calculated and ob- 
served images: 



Xred jy X/ 
0,y) 



[lZ hs {x,y)-r v {x,y)\ 



(27) 



where N is the number of degrees of freedom of the 
model, which is equal to the difference between the 
number of observational points and the number of 
model parameters, if the model is a linear function 
of these parameters. If the model is nonlinear, it is 
difficult (or even impossible) to determine the num- 
ber of degrees of freedom, since the computation 
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of one model becomes fairly resource-intensive [34j . 
We took N to be equal to the difference between 
the number of map pixels (60x60=3600 pixels for 
a 3"x3" map) and the number of free parame- 
ters (< 10). 

The reconstruction of the disk properties is re- 
duced to searching for the minimum X^cd as a mnc_ 
tion of the free model parameters. It is quite natural 
to single out four groups of parameters describing 
the following disk properties: 

I. Density distribution: 

• R ln — radius of the inner dust-depleted 
region; 

• -Rout — outer radius of the disk; 

• SQ as — gas surface density at a radius of 
1 AU; 

• p — index of the power-law function de- 
scribing the surface density. 

II. Thermal structure: 

• /istar — part of star radiative energy in- 
tercepted by the disk; 

• M star M — product of the stellar mass 
M s tar and the accretion rate M, describ- 
ing the accretion heating. 

III. Disk position: 

• i — inclination between the plane of the 
disk and the line of sight (i = 90° corre- 
sponds to an edge-on disk); 

• PA — position angle of the major axis of 
the disk image relative to the x axis of 
the map; 

• D — distance to the disk; 

• x s — shift of the central star position rel- 
ative to the map point (0,0) along the x 
axis; 

• y s — shift of the central star position rel- 
ative to the map point (0,0) along the y 
axis. 

IV. Dust particle properties: 

• — mass fraction of silicate grains mix- 
tured with graphite grains; 

• c m in — minimum size of the dust grains; 

• flmax — maximum size of the dust grains; 

• apow — index of the power-law dust grain 
size distribution. 

Some of these parameters can be kept fixed. Tests 
have shown that synthetic maps depend only weakly 
on the accretion rate. This is due to the fact that 



viscous heating dominates over heating by the stellar 
radiation in the inner, dense regions that are heated 
more, and whose maximum radiation occurs in the 
middle IR. The outer rarefied and cool regions of the 
disk are primarily responsible for the radiation in 
the millimeter. Therefore, it is not surprising that 
synthetic millimeter maps are not sensitive to M. 
Further, we can independently determine the posi- 
tion angle PA of the major axis of the disk image in 
the plane of the sky and the shifts of the image cen- 
ters x s , y s relative to the central star. We fixed the 
distance to the disk to be D = 140 pc, since there is 
reason to believe that the disk in CB 26 is a member 
of the Taurus- Auriga star- forming region [21]. It is 
reasonable to take the same parameters for the sili- 
cate and graphite particle-size distributions, since it 
has been shown that their evolutions do not differ 
appreciably [Jiirgen Blum, priv. comm.]. We are 
going to determine to what degree the dust parti- 
cles can grow, i.e., to determine their maximum size 
« max ; thus, for simplicity, we fixed the minimum size 
to be a m ; n = 5 • 10~ 7 cm and the distribution slope 
to be a pow = —3.5, corresponding to the interstellar 
medium (35j . 

After the above parameters are excluded from con- 
sideration, the following set of eight free parame- 
ters remains: (i? in , i? out , Sg as ,_p; fL stal ; i; f Si , a max ). 
If these parameters were determined using only spa- 
tially unresolved observations, this would yield cer- 
tain difficulties due to the small number of degrees 
of freedom in the model. The coverage of the fre- 
quency range considered would be determined by 
only about 20 points, corresponding to the available 
observations (HST, IRAS, Spitzer, Herschel, SMA, 
SCUBA); therefore, attempts to fit an SED alone 
often results in significant degeneracy of the model 
parameters. However, interferometric data provide 
additional information about the object's structure, 
making it possible to reduce the degeneracy, or, for 
a number of parameters, even remove this problem 
altogether. 

We used the conjugate-gradient method of Pow- 
ell [36[ to search for the minimum (l2"T]) .This is a 
fairly rapid algorithm that can solve the minimiza- 
tion problem efficiently if the function considered has 
"broad valleys". These regions are characterized by 
shallow, extended minima, such as are expected in 
the problem considered here. For example, the radi- 
ation flux from the disk is proportional to the disk 
mass and the Planck function of the dust temper- 
ature, if the disk is optically thin at a given wave- 
length. The same flux could come from both a cool, 
massive disk and a hotter, less massive disk. There- 
fore, the problem is degenerate for the disk mass and 
mean dust temperature in the disk, if the observa- 
tions are not spatially resolved. If interferometric 
observations are analyzed, this problem can result 
in a mutual dependence between the parameters for 
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Figure 1. An example of searching for the best-fit model using the Powell method. 



groups I and II. One of our goals was to study the 
degree of degeneracy of these parameters for spa- 
tially resolved observations at a level of 1" or better. 
The conjugate-gradient method is optimal for solv- 
ing such problems. Figure [1] presents an example of 
the convergence history of this method for one ob- 
served map. The computation time for one model is 
5-20 min using one 3 GHz processor. 

We determined the parameter uncertainties as fol- 
lows. The boundaries of a confidence interval for 
a given parameter were taken to be the points at 
which the value of Xrcd increases by unity compared 
to the minimum value, if all the other parameters are 
fixed and correspond to the minimum (37| . Thus, if 

min X?cd = Xrcd ( a i 0) > a 2 0) > a 3 0) >-)' then the C0Ilfi - 

dence interval for the parameter a\ is the interval 
[a^ — cr~~ + <7 + ^J , such that 

xL(ai 0) ±^4°\4V.) =min Xr 2 cd + l. (28) 

The non-linearity of the model (more exactly, the 
complexity of determining the number of degrees 
of freedom TV) means that these confidence inter- 
val cannot be taken to be accurate la uncertainties, 
but we consider them to be a measure of the param- 
eter uncertainty, close to the frequently used level of 
68.3%. 



6. RESULTS 

6.1. Best-Fit Models for Individual Maps 

In the first step of our study, we searched for best- 
fit models for each of the three maps independently. 
Figure [5] presents the results of this search, where 



the disk images are given at 110 GHz (left), 230 GHz 
(center), and 270 GHz (right). The upper row shows 
ideal images that would be observed with a tele- 
scope having a point-like beam. A prominent fea- 
ture is the presence of the central region devoid of 
dust, with radii of approximately 55 AU at 110 GHz 
and 35 AU at 230 and 270 GHz. The maximum ra- 
diation intensity is reached at the edge of this re- 
gion: 35, 280, and 560 mJy/arcsec 2 at 110, 230, 
and 270 GHz, respectively. The synthetic maps ob- 
tained by convolving the ideal images with a beam 
are presented in the second row. A comparison be- 
tween the ideal and convolved maps illustrates the 
difficulties arising in attempts to establish the disk 
sizes (morphology) and orientations using the corre- 
sponding observed maps. The maximum intensities 
in the synthetic maps are 10, 20, and 110 mJy/T^ 
at 110, 230, and 270 GHz, respectively. The inten- 
sities are expressed in terms of the effective beam 
solid angles f2„, for convenience in comparing with 
the observational data in the third row in Figure [5] 
(rJnQ = 1.63 arcsec 2 , Q230 = 0.16 arcsec 2 , Q270 = 
0.95 arcsec 2 ). The fourth row in Figure [5] presents 
the distribution of Xrcd (|27|) over the image. On the 
whole, the differences in the 230 and 270 GHz maps 
are distributed randomly, while the disk image seems 
to display structure at 110 GHz. This may provide 
evidence for the existence of systematic difference 
between the theoretical and observed maps at this 
frequency. 

Numerical values of some model parameters corre- 
sponding to the presented synthetic maps are given 
in Table[2] together with the confidence intervals cal- 
culated according to (|2"5)) and the minimum values 
Xrcd- Note the high level of agreement between the 
observed and synthetic maps: the mean deviation 
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Figure 2. Maps of the disk radiation intensity at 110 GHz (left column), 230 GHz (central column), and 270 GHz (right 
column). The upper row shows the ideal disk images for observations using a telescope with a point-like beam, the second 
row is the ideal disk images convolved with the beam, the third row is the observed images, and the fourth row is the maps 
of the difference between the observed (O) and model (C) images in noise units a v : ((O — C)/ov) 2 . 



over a map is less than 1.5 times the noise. 



We have analyzed in detail Xrod as a function 
of a number of important parameters: i?i n , i and 
a lllax . Figure [3] shows the dependence Xred(-^ia) m 
the vicinity of the minimum, keeping the remain- 



ing parameters fixed. In this case also, the param- 
eters obtained for 110 GHz differ appreciably from 
those obtained for the other two frequencies. The 
inner disk radius is approximately 37 ± 15 AU at 
230 and 270 GHz, and approximately 55 ± 5 AU at 



10 



Table 2. Parameters of the best-fit model. 



Parameter 


110 GHz 


230 GHz 


270 GHz 


230+270 GHz 


110+230+270 GHz 


Xred 

Rin, AU 

Rout, AU 
Eg", g/cm 2 
P 

i, deg 
/•£* 

/Si 

Qmax) Cm 
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54 + 3 
157+? 2 
782+^ 
-0 64+ 002 

79+12 
n n o+0.004 

u.uo_ 00g 
0.381°;°* 
< 1.7- 10" 2 
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78if 3 
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1.1 

36+^ 

Ju — 13 

210t^ 

7 r 7 +142 
'°' -126 

-0.87+°;°* 

78+ 9 21 
12+ 006 
47+° 18 

< 1.0 • 10~ 2 


1.7 

— 14 

222+g 

7-.Q+152 
' lu -147 
-0 81 + ' 04 

78l? 6 
0.081S;" 
0.59+°;« 
< 1.3 ■ 10~ 2 


2.2 

33+12 

lJJ — 15 

207 +2 J 

-0.69«;g 3 3 

82l? 7 
o oi<s+ - 005 
07+0.13 

u '°' -0.13 

< 1.7- 10" 2 
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Figure 3. Fit criterion X 2 ec j 38 a function of the 
inner disk radius Ri n . The other parameters were 

fixed at the values reached at the minimum 
Xq = mm X 2 c d- The arrows indicate the confidence 
intervals corresponding to X 2 e( j = Xo + b an< l the 
points on the arrows correspond to the minimum. 



110 GHz. The presence of a dust-free region with a 
size of 45 ± 5 AU was suggested in [24| based on the 
presence of a central plateau in the 230 GHz map. 
Figure [3] demonstrates that neither the new obser- 
vations obtained u sing the PdBI (110 and 230 GHz) 
nor those used in [2J] (obtained using the SMA at 
270 GHz) can be explained without assuming the 
existence of a large, free-dust central region. This 
region could have arisen due to both dynamical ef- 
fects, such as the sweeping up of dust by a binary 
star, and evolutionary effects, such as the formation 
of planetesimals in the central regions. In this case, 
the absence of radiation from the central region can 
be explained by inefficient reprocessing of the radi- 
ation from the central star, rather than an absence 
of matter. Such holes are often observed in other 
disks [H 13. 

The disk inclination i is an important parameter 
for modeling the bipolar outflow from CB 26. The 
accuracy in this parameter can be estimated using 
Figure [H which presents the dependence X?ed(*) f° r 
the three frequencies. We show confidence intervals 
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Figure 4. Fit criterion X 2 e d ^ a function of the 
inclination angle i of the disk to the plane of the sky. 

The other parameters were fixed at the values 
reached at the minimum Xo = mm X 2 c( j- The arrows 
indicate the confidence intervals corresponding to 
X 2 cd = Xo + 1 an d * < 90° (this orientation follows 
from observations in molecular lines). The points on 
the arrows correspond to the minimum. For i < 90° , 
the south pole of the disk axis is oriented toward the 
Earth. 

for i < 90°, since the southern edge of the disk is 
nearer to the Earth, as is suggested by observations 
of the outflow in molecular lines. The mean incli- 
nation over the three maps is 76°. None of these 
millimeter maps are sensitive to the edge of the disk 
that is nearest the observer (the symmetry of X^cd(*) 
relative to 90°). 

Figure [6] shows the dependence x;? ed (a max ). The 
synthetic maps are not sensitive to the maximum 
size of the dust particles in the range from interstel- 
lar values to 10 -2 cm. Our model suggests an upper 
limit to the maximum size of the dust particles in 
CB 26 of 0.02 cm. 

6.2. Combined Model 

The fifth column in Table [2] presents the model pa- 
rameters determined using the 230 and 270 GHz 
maps jointly. The maps at these two frequencies ob- 
tained with different angular resolutions using differ- 
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Figure 5. Dust temperature (left) and gas density (right) distributions for the disk models whose best-fit parameters 

correspond to the observed maps at 230 and 270 GHz. 



110 GHz 

230 GHz 

270 GHz 




1e-05 0.0001 0.001 0.01 

a max [cm] 



*-rcd 



maximum size of the dust particles a max - The other 
parameters were fixed at the values reached at the 
minimum \q = mm Xr C d' arrows indicate the 

confidence intervals corresponding to Xr C d = Xo + 1- 



ent instruments yield quite similar disk parameters, 
whereas the 110 GHz map gives different parame- 
ters and shows the presence of correlated features in 
the residuals (Figure [5]). Therefore, we excluded the 
110 GHz map from our final search for the disk pa- 
rameters. The mean value of Xrcd = (X230 + X27o)/2 £n as and the portion of the stellar radiation inter- 



/L* = 0.08L Q . The full optical depth of the disk 
is close to unity at all three frequencies. The maxi- 
mum optical depth is reached in the vicinity of the 
disk inner boundary, and is equal to 0.35, 0.6, and 
0.9 for 110, 230, and 270 GHz, correspondingly. 



6.3. Model Degeneracy 



The fact that we observe the disk edge-on and can- 
not see the direct radiation from the star hinders de- 
termination of the spectral type of the central star. 
The lack of information about the effective temper- 
ature of the star means that the thermal structure 
of the disk must be reconstructed, and the charac- 
teristic dust temperature required to determine the 
disk mass [l(| must be estimated. There is a fun- 
damental (physical) constraint on the joint, simul- 
taneous determination of the mass and characteris- 
tic temperature of the disk in the case of spatially 
unresolved observations in the millimeter. We con- 
structed maps of Xred as functions of the structural 
parameters in order to study the degeneracy of spa- 
tially resolved observations. Figure shows x re d as 
functions of the normalization of the surface density 



v-red 

was used for x re d- The resUi ts obtained using all 
three maps are given in the last column of Table [2] 
for comparison. 

The physical structure of the disk corresponding 
to the joint minimum Xred ^ s snown m Figure|5] The 
dust temperature varies from 10 K in the peripheral 
regions near the equatorial plane to 100 K in the at- 
mosphere of the disk at its center. The maximum gas 
density in the equatorial plane reaches 10~ 12 g/cm 3 . 
The normalized surface density, power-law index, 
and portion of the stellar radiation intercepted by 
the disk are Sf s = 710 g/cm 2 , p = -0.8, and 



cepted by the disk /L*. For ease of viewing, the 
latter quantity has been transformed to the effective 
temperature of the star using the standard formula 
fL* = f^R 2 star ■ a B T s 4 tal . for R stal = R J = 0.1. 
The values of parameters from the dark-grey re- 
gion with Xicd ^ 2 (Figure [7]) describe the obser- 
vations equally well. This region of parameters can 
be described by the law Ep as T Q = const, where 
a = 1.3-r- 1.7. x^d behaves similarly in the (T,Q &s ,p) 
parameter plane (Figure |8]), with the steeper de- 
crease in the surface density toward the periphery 
corresponding to greater values of S| as . 
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Figure 7. Fit criterion X^ed 33 a f unc ti° n of the 
normalization of the surface density SQ as and the 
temperature T s t a r of the central star. 



Figure 8. Fit criterion Xr C d as a function of the 
normalization of the surface density SQ as and the 
index p of the power-law function for Ex- 



6.4. Origin of the Differences of 110 GHz Map 

All the available data indicate that the 110 GHz map 
(corresponding to longest of the wavelengths con- 
sidered) differs from the remaining two. First, the 
residuals display correlated features within the disk 
image (Figure^, suggesting a systematic difference 
between the observational and theoretical data at 
this frequency, or that the derived Xrcd represents 
a local minimum. Second, the disk parameters de- 
rived using this map differ appreciably from those 



obtained using the 230 and 270 GHz maps (Tabled . 
Third, if the point corresponding to 110 GHz is plot- 
ted on the SED (Figure it deviates appreciably 
from the straight line in the Rayleigh- Jeans part of 
the s pec trum. To plot Figure Owe used the fluxes 
from [24| together with data from the Herschel. The 
deviation of a linear SED and the differences in the 
disk parameters at 110 GHz and at the other fre- 
quencies can be explained in the following way: 

• there exists an additional radiative mechanism 
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Figure 9. SED of CB 26. The open circles show the data from which the envelope's contribution has been 
subtracted. For ease of viewing, the straight line shows the slope of the Rayleigh-Jeans part of the spectrum, derived 

using the last four points, except for 110 GHz. 



that is not taken into account in the model, 
which is not important at high frequencies, but 
becomes appreciable at 110 GHz (e.g., free- 
free radiation arising in the accretion region 
or bipolar outflow); 

• an envelope around the disk could also con- 
tribute to the radiation, since the 110 GHz 
beam is larger than the 230 and 270 GHz 
beams (we are planning to use CO line ob- 
servations to estimate the envelope's contribu- 
tion); 

• it is difficult to establish an accurate dust 
model (chemical composition and extinction 
efficiency coefficients). 

These (or other) effects could result in the overall 
difference between the 110 GHz map and the maps 
at 230 and 270 GHz. It seems reasonable to use the 
results based on the latter two maps (rather than all 
three) as a final disk model, until we determine the 
nature of the operating additional mechanism and 
take it into account. 



7. CONCLUSIONS 

The appearance of modern radio interferometers 
with high sensitivities and good angular resolutions 
(SMA, Plateau de Bure, CARMA, ALMA) has made 



it possible to obtain spatially resolved images of pro- 
toplanetary disks. The reconstruction of the phys- 
ical structures of the disks from these observations 
is a challenging task that requires an integrated ap- 
proach, that is simultaneously self-consistent physi- 
cally and correct mathematically, aimed at modeling 
and fitting the observations. One of the main goals 
of this paper was to develop such an approach. 

The key characteristic of the protoplanetary disk 
model we used to reconstruct the disk parameters 
is a compromise between the completeness of the 
description provided and the time required to com- 
pute the physical structure of the disk. The use of 
a two-frequency approximation and the correspond- 
ing mean opacities in the radiative transfer model, 
as well as the iterative scheme employed to solve 
the equation of hydrostatic equilibrium, yield a short 
model computation time. Together with the rapid 
computation of the theoretical maps, the adaptive 
algorithm for convolution of the images with the 
beam, and the efficient method used to minimize 
Xred' this made it possible to use the developed code 
to search for best-fit model parameters of the disk. 

We applied this code to determine the physical 
structure of the protoplanetary disk in CB 26 us- 
ing interferometric images in the millimeter. Obser- 
vational data at 110, 230, and 270 GHz from the 
SMA, IRAM Plateau de Bure, and OVRO inter- 
ferometers were used to search for best-fit models. 
Best-fit model parameters were obtained for each of 
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the frequencies, for all the frequencies jointly, and 
for 230 and 270 GHz jointly. All three observational 
maps suggest the existence of a central, dust-free 
region in the disk, approximately 35 AU in radius. 
This corresponds well to the value of 45 AU obtained 
in [24| at 230 GHz. However, our simulated maps 
yield a disk inclination of 78°, which is lower than 
the value of 85° used in (23j. This should lead to 
an appreciable decrease in the bipolar outflow ve- 
locity for the model [23[, and to a lower rate of an- 
gular momentum transfer from the disk. We cannot 
determine whether or not the dust particles in the 
disk can grow appreciably in size (see also (24||); we 
can only establish an upper bound on the dust size 
(omax < 0.02 cm) using our dust model, comprising 
a mixture of silicate and graphite particles with a 
power-law size distribution. 

We came across the following difficulties when 
simulating CB 26. The 110 GHz map stands out 
from the others due to the change in the SED slope, 
the presence of symmetrical features in the residual 
map within the disk image, and differences in the 
physical parameters implied by the best-fit model. 
These may indicate that (a) the model does not in- 
clude free-free radiation, although this may be com- 
parable with the thermal radiation of the disk at 
longer wavelengths; (b) the contribution from a disk 
envelope has been neglected; (c) the dust model is 
not sufficiently realistic. All these factors must be 
considered in more detail. Our analysis also indi- 



cates a wide region of degeneracy between the ther- 
mal and density characteristics of the disk. In our 
opinion, this is a fundamental problem hindering our 
ability to unambiguously establish the disk param- 
eters. This problem can be avoided in objects, in 
which the radiation of the central star is not screened 
by the circumstellar disk, and the stellar parameters 
can be determined independently. It is also obvious 
that an increase in angular resolution is necessary for 
this degeneracy to be removed. We plan to study in 
future whether or not the region of parameter de- 
generacy can be diminished using maps from other 
wavebands. 
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